#Analysis of new mapping 15/3/11
#maps all populations simultaneously and then splits reads

#see commands_2011_03_09.txt
#Read in AllDiffs files

grouse1mar <- read.table('/Paterson/Datafiles/grouse/newbler_mapping/grouse1_all_diffs_2011_03_14.txt',header= TRUE, stringsAsFactors=FALSE,strip.white=TRUE,sep="\t")
grouse2mar <- read.table('/Paterson/Datafiles/grouse/newbler_mapping/grouse2_all_diffs_2011_03_14.txt',header= TRUE, stringsAsFactors=FALSE,strip.white=TRUE,sep="\t")
grouse3mar <- read.table('/Paterson/Datafiles/grouse/newbler_mapping/grouse3_all_diffs_2011_03_14.txt',header= TRUE, stringsAsFactors=FALSE,strip.white=TRUE,sep="\t")

grouse.mar <- grouse1mar
names(grouse.mar)[names(grouse.mar) %in% c('ref.reads','var.reads')] <- c('pop1.ref.reads','pop1.var.reads')
grouse.mar$pop1.depth <- rowSums(grouse.mar[,c('pop1.ref.reads','pop1.var.reads')])

grouse.mar <- cbind(grouse.mar,grouse2mar[,c('ref.reads','var.reads')])
names(grouse.mar)[names(grouse.mar) %in% c('ref.reads','var.reads')] <- c('pop2.ref.reads','pop2.var.reads')
grouse.mar$pop2.depth <- rowSums(grouse.mar[,c('pop2.ref.reads','pop2.var.reads')])

grouse.mar <- cbind(grouse.mar,grouse3mar[,c('ref.reads','var.reads')])
names(grouse.mar)[names(grouse.mar) %in% c('ref.reads','var.reads')] <- c('pop3.ref.reads','pop3.var.reads')
grouse.mar$pop3.depth <- rowSums(grouse.mar[,c('pop3.ref.reads','pop3.var.reads')])

grouse.mar$pop12.ref.reads <- rowSums(grouse.mar[,c('pop1.ref.reads','pop2.ref.reads')])
grouse.mar$pop12.var.reads <- rowSums(grouse.mar[,c('pop1.var.reads','pop2.var.reads')])
grouse.mar$pop12.depth <- rowSums(grouse.mar[,c('pop12.ref.reads','pop12.var.reads')])

#need to produce something that looks like grouse.df3
grabSNPs <- function(contig.map.df, grouse.mar,cDNA){
	require(GeneR)
	#contig.map.df is one of the datafiles held in contig.map
	contig.map.df$contig <- sub('\\(\\w*)','\1',contig.map.df[,1])
	out.df <- data.frame(cDNA=character(0),id=character(0),contig=character(0),
		super.start.pos=numeric(0),start.pos=numeric(0),end.pos=numeric(0),
		ref.nucl=character(0),var.nucl=character(0),
		pop1.ref.reads=numeric(0), pop1.var.reads=numeric(0),pop1.depth=numeric(0),
		pop2.ref.reads=numeric(0), pop2.var.reads=numeric(0),pop2.depth=numeric(0),
		pop2.ref.reads=numeric(0), pop3.var.reads=numeric(0),pop3.depth=numeric(0),
		pop12.ref.reads=numeric(0), pop12.var.reads=numeric(0),pop12.depth=numeric(0))
	
		
	for(i in 1:nrow(contig.map.df)){
		contig.name <- contig.map.df$contig[i]
		tmp.mar <- grouse.mar[grouse.mar$Reference==contig.name,]
		if(nrow(tmp.mar)==0) next
		tmp.mar$cDNA <- cDNA
		tmp.mar$contig <- contig.name
		if(contig.map.df$direction[i]=="forward"){
			tmp.mar$super.start.pos <- contig.map.df[i,'concat.start'] + tmp.mar$start.pos - 1
			}else{
			tmp.mar$super.start.pos <- contig.map.df[i,'concat.end'] - tmp.mar$end.pos - 99
			tmp.mar$ref.nucl[tmp.mar$ref.nucl!="-"] <- sapply(tmp.mar$ref.nucl[tmp.mar$ref.nucl!="-"],strComp)
			tmp.mar$var.nucl[tmp.mar$var.nucl!="-"] <- sapply(tmp.mar$var.nucl[tmp.mar$var.nucl!="-"],strComp)
			}
		tmp.mar$id <- paste(cDNA,tmp.mar$super.start.pos,contig.name,tmp.mar$start.pos,tmp.mar$ref.nucl,tmp.mar$var.nucl,sep="_")		
		out.df <- rbind(out.df, tmp.mar[,c('cDNA','id','contig','super.start.pos','start.pos','end.pos',
			'ref.nucl','var.nucl','pop1.ref.reads','pop1.var.reads','pop1.depth',
			'pop2.ref.reads','pop2.var.reads','pop2.depth','pop3.ref.reads','pop3.var.reads','pop3.depth',
			'pop12.ref.reads','pop12.var.reads','pop12.depth')])

		}
	out.df
	}
#debug(grabSNPs)
#tst <- grabSNPs(contig.map[['cDNA_260-1']],grouse.mar,'cDNA_260-1')

grouse.mar.df <- data.frame(cDNA=character(0),id=character(0),contig=character(0),
		super.start.pos=numeric(0),start.pos=numeric(0),end.pos=numeric(0),
		ref.nucl=character(0),var.nucl=character(0),
		pop1.ref.reads=numeric(0), pop1.var.reads=numeric(0),pop1.depth=numeric(0),
		pop2.ref.reads=numeric(0), pop2.var.reads=numeric(0),pop2.depth=numeric(0),
		pop2.ref.reads=numeric(0), pop3.var.reads=numeric(0),pop3.depth=numeric(0),
		pop12.ref.reads=numeric(0), pop12.var.reads=numeric(0),pop12.depth=numeric(0))

for(cdna in names(contig.map)){
	if(identical(contig.map[[cdna]],0)) {
		cat(cdna,' skipped\t')
		next
		}
	grouse.mar.df <- rbind(grouse.mar.df,grabSNPs(contig.map[[cdna]],grouse.mar,cdna))
	cat(cdna,' read\t')
	}
#41457 rows as output

grouse.mar.df$pop123.ref.reads <- rowSums(grouse.mar.df[,c('pop12.ref.reads','pop3.ref.reads')])
grouse.mar.df$pop123.var.reads <- rowSums(grouse.mar.df[,c('pop12.var.reads','pop3.var.reads')])
grouse.mar.df$pop123.depth <- rowSums(grouse.mar.df[,c('pop12.depth','pop3.depth')])
grouse.mar.df$ref.nucl <- as.character(grouse.mar.df$ref.nucl)
grouse.mar.df$var.nucl <- as.character(unlist(grouse.mar.df$var.nucl))
#Run a filter for good quality snps
#total depth >30 and var count > 5

grouse.mar.df2 <- grouse.mar.df[grouse.mar.df$pop123.depth > 30 & grouse.mar.df$pop123.var.reads >5,]
#19974
grouse.mar.df2 <- grouse.mar.df2[apply(X=grouse.mar.df2[,c('pop1.depth','pop1.depth','pop1.depth')],1,function(X){all(X>0)}),]
#19974

#trim out 2 snps that appear at the same start position
tmp.filt <- data.frame(snp.start=paste(grouse.mar.df2$contig,grouse.mar.df2$start.pos,sep="_"),stringsAsFactors=FALSE)
tmp.filt$duplicate <- tmp.filt$snp.start %in% names(table(tmp.filt$snp.start))[table(tmp.filt$snp.start)>1]

grouse.mar.df2 <- grouse.mar.df2[!tmp.filt$duplicate,]
#18865

grouse.mar.df2$p.hat.all <- apply(X=grouse.mar.df2[,c('pop1.var.reads','pop2.var.reads','pop3.var.reads', 'pop1.depth','pop2.depth','pop3.depth')],1,function(X){
	tmp.X <- X[c('pop1.var.reads','pop2.var.reads','pop3.var.reads')]/X[c('pop1.depth','pop2.depth','pop3.depth')]
	mean(tmp.X)
	})
grouse.mar.df2$p.hat.12v3 <- apply(X=grouse.mar.df2[,c('pop12.var.reads','pop3.var.reads', 'pop12.depth','pop3.depth')],1,function(X){
	tmp.X <- X[c('pop12.var.reads','pop3.var.reads')]/X[c('pop12.depth','pop3.depth')]
	mean(tmp.X)
	})
grouse.mar.df2$p.hat.1v2 <- apply(X=grouse.mar.df2[,c('pop1.var.reads','pop2.var.reads', 'pop1.depth','pop2.depth')],1,function(X){
	tmp.X <- X[c('pop1.var.reads','pop2.var.reads')]/X[c('pop1.depth','pop2.depth')]
	mean(tmp.X)
	})
grouse.mar.df2$p.hat.12v3 <- apply(X=grouse.mar.df2[,c('pop12.var.reads','pop3.var.reads', 'pop12.depth','pop3.depth')],1,function(X){
	tmp.X <- X[c('pop12.var.reads','pop3.var.reads')]/X[c('pop12.depth','pop3.depth')]
	mean(tmp.X)
	})


########################################
## Calculate FST

ifst <- seq(from=0.001,to=0.99,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df2[grouse.mar.df2$p.hat.all >0 &grouse.mar.df2$p.hat.all < 1,]),ncol=100,dimnames=list(grouse.mar.df2$id[grouse.mar.df2$p.hat.all >0 &grouse.mar.df2$p.hat.all < 1],paste("theta",itheta,sep="_")))


#grouse.mar.df2[1,]
#LikBetaBinomialPops(depth=grouse.mar.df2[1,c('pop1.depth','pop1.depth','pop1.depth')],
#	count=grouse.mar.df2[1,c('pop1.var.reads','pop1.var.reads','pop1.var.reads')],
#	p.hat=grouse.mar.df2[1,'p.hat.all'], theta= 0.7)

lik.theta.list <- vector(mode='list',length=3)
names(lik.theta.list) <- c('all.pops','pop1v2','pop12v3')

for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df2[grouse.mar.df2$p.hat.all >0 &grouse.mar.df2$p.hat.all < 1,c('pop1.var.reads','pop2.var.reads','pop3.var.reads', 'pop1.depth','pop2.depth','pop3.depth','p.hat.all')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop1.depth','pop2.depth','pop3.depth')],
		count=X[c('pop1.var.reads','pop2.var.reads','pop3.var.reads')],
		p.hat= X['p.hat.all'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 0.071
ifst[null.which]

#cut and paste from previous analysis: beta-binomial method
require(multtest)
BBlikTest.mar.df <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(lik.theta.fit)))
BBlikTest.mar.df$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df$chi.sq <- 2*(BBlikTest.mar.df$max.lik-BBlikTest.mar.df$null.lik)
BBlikTest.mar.df$p.val <- 1-pchisq(BBlikTest.mar.df$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df$p.val,'Bonferroni')
BBlikTest.mar.df$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df$p.val,'Hochberg')
BBlikTest.mar.df$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled <- aggregate(lik.theta.fit,by=list(grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.all >0 &grouse.mar.df2$p.hat.all < 1]),sum)
BBlikTest.mar.pooled.df <- data.frame(cDNA=BBlikTest.mar.pooled[,1],null.lik = BBlikTest.mar.pooled[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled)))
BBlikTest.mar.pooled.df$cDNA <- as.character(BBlikTest.mar.pooled.df$cDNA)	
BBlikTest.mar.pooled.df$max.lik <- apply(BBlikTest.mar.pooled[,-1],1,max)
BBlikTest.mar.pooled.df$chi.sq <- 2*(BBlikTest.mar.pooled.df$max.lik-BBlikTest.mar.pooled.df$null.lik)
BBlikTest.mar.pooled.df$p.val <- 1-pchisq(BBlikTest.mar.pooled.df$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df$p.val,'Bonferroni')
BBlikTest.mar.pooled.df$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df$p.val,'Hochberg')
BBlikTest.mar.pooled.df$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df$fst <- 0

BBlikTest.mar.pooled.df$fst <- apply(BBlikTest.mar.pooled[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.99,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df$fst <- as.numeric(BBlikTest.mar.pooled.df$fst)
BBlikTest.mar.pooled.df <- BBlikTest.mar.pooled.df[!is.na(BBlikTest.mar.pooled.df$null.lik),]

#some big differences with population 3 seen, essentially private alleles

lik.theta.list[['all.pops']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df,BBlikTest.pooled.df=BBlikTest.mar.pooled.df,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.all >0 &grouse.mar.df2$p.hat.all < 1])


#try to repeat analysis just on pop1 vs pop2

ifst <- seq(from=0.001,to=0.5,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df2[grouse.mar.df2$p.hat.1v2 >0 &grouse.mar.df2$p.hat.1v2 < 1,]),ncol=100,dimnames=list(grouse.mar.df2$id[grouse.mar.df2$p.hat.1v2 >0 &grouse.mar.df2$p.hat.1v2 < 1],paste("theta",itheta,sep="_")))


#grouse.mar.df2[1,]
#LikBetaBinomialPops(depth=grouse.mar.df2[1,c('pop1.depth','pop1.depth','pop1.depth')],
#	count=grouse.mar.df2[1,c('pop1.var.reads','pop1.var.reads','pop1.var.reads')],
#	p.hat=grouse.mar.df2[1,'p.hat.all'], theta= 0.7)

for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df2[grouse.mar.df2$p.hat.1v2 >0 &grouse.mar.df2$p.hat.1v2 < 1,c('pop1.var.reads','pop2.var.reads', 'pop1.depth','pop2.depth','p.hat.1v2')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop1.depth','pop2.depth')],
		count=X[c('pop1.var.reads','pop2.var.reads')],
		p.hat= X['p.hat.1v2'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 
ifst[null.which] # 0.011


BBlikTest.mar.df1v2 <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(lik.theta.fit)))
BBlikTest.mar.df1v2$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df1v2$chi.sq <- 2*(BBlikTest.mar.df1v2$max.lik-BBlikTest.mar.df1v2$null.lik)
BBlikTest.mar.df1v2$p.val <- 1-pchisq(BBlikTest.mar.df1v2$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df1v2$p.val,'Bonferroni')
BBlikTest.mar.df1v2$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df1v2$p.val,'Hochberg')
BBlikTest.mar.df1v2$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled.1v2 <- aggregate(lik.theta.fit,by=list(grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.1v2 >0 &grouse.mar.df2$p.hat.1v2 < 1]),sum)
BBlikTest.mar.pooled.df1v2 <- data.frame(cDNA=BBlikTest.mar.pooled.1v2[,1],null.lik = BBlikTest.mar.pooled.1v2[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled.1v2)))
BBlikTest.mar.pooled.df1v2$cDNA <- as.character(BBlikTest.mar.pooled.df1v2$cDNA)	
BBlikTest.mar.pooled.df1v2$max.lik <- apply(BBlikTest.mar.pooled.1v2[,-1],1,max)
BBlikTest.mar.pooled.df1v2$chi.sq <- 2*(BBlikTest.mar.pooled.df1v2$max.lik-BBlikTest.mar.pooled.df1v2$null.lik)
BBlikTest.mar.pooled.df1v2$p.val <- 1-pchisq(BBlikTest.mar.pooled.df1v2$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df1v2$p.val,'Bonferroni')
BBlikTest.mar.pooled.df1v2$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df1v2$p.val,'Hochberg')
BBlikTest.mar.pooled.df1v2$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df1v2$fst <- 0

BBlikTest.mar.pooled.df1v2$fst <- apply(BBlikTest.mar.pooled.1v2[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.5,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df1v2$fst <- as.numeric(BBlikTest.mar.pooled.df1v2$fst)
BBlikTest.mar.pooled.df1v2 <- BBlikTest.mar.pooled.df1v2[!is.na(BBlikTest.mar.pooled.df1v2$null.lik),]

lik.theta.list[['pop1v2']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df1v2,BBlikTest.pooled.df=BBlikTest.mar.pooled.df1v2,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.1v2 >0 &grouse.mar.df2$p.hat.1v2 < 1])


#should also contrast pop12 vs pop3

ifst <- seq(from=0.001,to=0.99,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df2[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1,]),ncol=100,dimnames=list(grouse.mar.df2$id[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1],paste("theta",itheta,sep="_")))


for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df2[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1,c('pop12.var.reads','pop3.var.reads', 'pop12.depth','pop3.depth','p.hat.12v3')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop12.depth','pop3.depth')],
		count=X[c('pop12.var.reads','pop3.var.reads')],
		p.hat= X['p.hat.12v3'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 
ifst[null.which] # 0.0609394


BBlikTest.mar.df12v3 <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(grouse.mar.df2[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1,])))
BBlikTest.mar.df12v3$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df12v3$chi.sq <- 2*(BBlikTest.mar.df12v3$max.lik-BBlikTest.mar.df12v3$null.lik)
BBlikTest.mar.df12v3$p.val <- 1-pchisq(BBlikTest.mar.df12v3$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df12v3$p.val,'Bonferroni')
BBlikTest.mar.df12v3$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df12v3$p.val,'Hochberg')
BBlikTest.mar.df12v3$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled.12v3 <- aggregate(lik.theta.fit,by=list(grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1]),sum)
BBlikTest.mar.pooled.df12v3 <- data.frame(cDNA=BBlikTest.mar.pooled.12v3[,1],null.lik = BBlikTest.mar.pooled.12v3[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled.12v3)))
BBlikTest.mar.pooled.df12v3$cDNA <- as.character(BBlikTest.mar.pooled.df12v3$cDNA)	
BBlikTest.mar.pooled.df12v3$max.lik <- apply(BBlikTest.mar.pooled.12v3[,-1],1,max)
BBlikTest.mar.pooled.df12v3$chi.sq <- 2*(BBlikTest.mar.pooled.df12v3$max.lik-BBlikTest.mar.pooled.df12v3$null.lik)
BBlikTest.mar.pooled.df12v3$p.val <- 1-pchisq(BBlikTest.mar.pooled.df12v3$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df12v3$p.val,'Bonferroni')
BBlikTest.mar.pooled.df12v3$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df12v3$p.val,'Hochberg')
BBlikTest.mar.pooled.df12v3$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df12v3$fst <- 0

BBlikTest.mar.pooled.df12v3$fst <- apply(BBlikTest.mar.pooled.12v3[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.99,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df12v3$fst <- as.numeric(BBlikTest.mar.pooled.df12v3$fst)
BBlikTest.mar.pooled.df12v3 <- BBlikTest.mar.pooled.df12v3[!is.na(BBlikTest.mar.pooled.df12v3$null.lik),]


lik.theta.list[['pop12v3']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df12v3,BBlikTest.pooled.df=BBlikTest.mar.pooled.df12v3,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df2$cDNA[grouse.mar.df2$p.hat.12v3 >0 &grouse.mar.df2$p.hat.12v3 < 1])


#need to perform bootstrapping, to permute the lik.theta.fit matrix vs cDNA for the aggregate function
#(move that over to the CGR server)

save(lik.theta.list,file='lik_theta_list.txt',ascii=TRUE)

#setup matrix to contain data
mat.pooled <- matrix(0,ncol=10,nrow=nrow(lik.theta.list[[1]]$BBlikTest.pooled.df))

for(i in 1:10){
	  
	  tmp.lik <- tmp.lik[sample(1:nrow(tmp.lik)),]
	tmp.agg <- aggregate(tmp.lik,by=list(lik.theta.list[[1]]$cDNA),sum)
	
	tmp.pooled <- data.frame(cDNA=tmp.agg[,1],null.lik = tmp.agg[,lik.theta.list[[1]]$null.which+1],
	max.lik = rep(0,nrow(tmp.agg)))
	tmp.pooled$cDNA <- as.character(tmp.pooled$cDNA)
	
	tmp.pooled$max.lik <- apply(tmp.agg[,-1],1,max)
	tmp.pooled$chi.sq <- 2*(tmp.pooled$max.lik-tmp.pooled$null.lik)
	tmp.pooled$p.val <- 1-pchisq(tmp.pooled$chi.sq,1)
	mat.pooled[,i] <- tmp.pooled$chi.sq
	}

tst.out <- data.frame(cDNA=tmp.pooled[,1],stringsAsFactors=FALSE)
tst.out <- cbind(tst.out,t(apply(X=mat.pooled,1,function(X){quantile(X,c(0.9,0.95,0.99,0.999))})))
save(tst.out,file='tst_out.txt',ascii=TRUE)

#see files in euler01:~/Datfiles/grouse/newbler_mapping/R_analysis

load("pop_all_out.txt") #pop.all
load("pop_1v2_out.txt")	#pop.1v2
load("pop_12v3_out.txt") #pop.12v3

#added extra trimming

#gets rid of indels
grouse.mar.df3 <- grouse.mar.df2[!apply(grouse.mar.df2[,c('ref.nucl','var.nucl')],1,function(X){any(X=="-")}),]
source("fst_mapping_code_df3.R")
save(lik.theta.list2,file='lik_theta_list2.txt',ascii=TRUE)

load("pop_all_out2.txt") #pop2.all
load("pop_1v2_out2.txt")	#pop2.1v2
load("pop_12v3_out2.txt") #pop2.12v3

#a bunch of code lost after laptop crashed
#used to output to excel files
#involved a bunch of merge commands to get bootstrapped p values

summary.fst <- read.table('../Stu_output/mar_output/summary_fst_2011_03_21.txt',header=TRUE,stringsAsFactors=FALSE,sep="\t")

############################
## Tajima D
############################

#sliding window 1000, min depth 50

grouse.TjD.list.mar <- vector('list',length=length(unique(grouse.mar.df3$cDNA)))
names(grouse.TjD.list.mar) <- unique(grouse.mar.df3$cDNA)

for(gene in unique(grouse.mar.df3$cDNA)){ #unique(grouse.df3$cDNA)
	tst.in <- grouse.mar.df3[grouse.mar.df3$cDNA==gene,]
	grouse.TjD.list.mar[[gene]][[1]] <- try(
		TajimaD(tst.in,sliding.window=1000,pos="super.start.pos",count="pop1.var.reads",depth="pop1.depth",n=50)
		)
	grouse.TjD.list.mar[[gene]][[2]] <- try(
		TajimaD(tst.in,sliding.window=1000,pos="super.start.pos",count="pop2.var.reads",depth="pop2.depth",n=50)
		)
	grouse.TjD.list.mar[[gene]][[3]] <- try(
		TajimaD(tst.in,sliding.window=1000,pos="super.start.pos",count="pop3.var.reads",depth="pop3.depth",n=50)
		)
	grouse.TjD.list.mar[[gene]][[4]] <- try(
		TajimaD(tst.in,sliding.window=1000,pos="super.start.pos",count="pop12.var.reads",depth="pop12.depth",n=50)
		)

	names(grouse.TjD.list.mar[[gene]]) <- c('pop1','pop2','pop3','pop12')
	}

grouse.TjD.mar.df <- data.frame(cDNA = unique(grouse.mar.df3$cDNA),stringsAsFactors=FALSE)
grouse.TjD.mar.df[,c('pop1.min','pop1.max','pop2.min','pop2.max','pop3.min','pop3.max','pop12.min','pop12.max')] <- 0

for(gene in 1:nrow(grouse.TjD.mar.df)){
	if(!is.null(grouse.TjD.list.mar[[gene]][[1]])&nrow(grouse.TjD.list.mar[[gene]][[1]])>0) grouse.TjD.mar.df[gene,c('pop1.min','pop1.max')] <- range(grouse.TjD.list.mar[[gene]][[1]]$D)
	if(!is.null(grouse.TjD.list.mar[[gene]][[2]])&nrow(grouse.TjD.list.mar[[gene]][[2]])>0) grouse.TjD.mar.df[gene,c('pop2.min','pop2.max')] <- range(grouse.TjD.list.mar[[gene]][[2]]$D)
	if(!is.null(grouse.TjD.list.mar[[gene]][[3]])&nrow(grouse.TjD.list.mar[[gene]][[3]])>0) grouse.TjD.mar.df[gene,c('pop3.min','pop3.max')] <- range(grouse.TjD.list.mar[[gene]][[3]]$D)
	if(!is.null(grouse.TjD.list.mar[[gene]][[4]])&nrow(grouse.TjD.list.mar[[gene]][[4]])>0) grouse.TjD.mar.df[gene,c('pop12.min','pop12.max')] <- range(grouse.TjD.list.mar[[gene]][[4]]$D)

	}

# Need to produce a randomisation test for Tajima D

grouse.TjD.mar.unlist <- data.frame(cDNA=character(0),pos=numeric(0),D1=numeric(0),S1=numeric(0), D2=numeric(0),S2=numeric(0),D3=numeric(0),S3=numeric(0),D12=numeric(0),S12=numeric(0))

for(gene in names(grouse.TjD.list.mar)){
	tmp.pos <- unique(c(grouse.TjD.list.mar[[gene]]$pop1$pos,
		grouse.TjD.list.mar[[gene]]$pop2$pos, grouse.TjD.list.mar[[gene]]$pop3$pos,
		grouse.TjD.list.mar[[gene]]$pop12$pos))
	if(length(tmp.pos)==0) next #jump out
	tmp.df <- data.frame(cDNA=gene,pos=tmp.pos,stringsAsFactors=FALSE,
	 	D1=as.numeric(rep(NA,length(tmp.pos))),S1=as.numeric(rep(NA,length(tmp.pos))),
	 	D2=as.numeric(rep(NA,length(tmp.pos))),S2=as.numeric(rep(NA,length(tmp.pos))),
	 	D3=as.numeric(rep(NA,length(tmp.pos))),S3=as.numeric(rep(NA,length(tmp.pos))),
	 	D12=as.numeric(rep(NA,length(tmp.pos))),S12=as.numeric(rep(NA,length(tmp.pos))))

	tmp.df[match(grouse.TjD.list.mar[[gene]]$pop1$pos,tmp.pos),c('D1','S1')] <- grouse.TjD.list.mar[[gene]]$pop1[,c('D','S')]
	tmp.df[match(grouse.TjD.list.mar[[gene]]$pop2$pos,tmp.pos),c('D2','S2')] <- grouse.TjD.list.mar[[gene]]$pop2[,c('D','S')]
	tmp.df[match(grouse.TjD.list.mar[[gene]]$pop3$pos,tmp.pos),c('D3','S3')] <- grouse.TjD.list.mar[[gene]]$pop3[,c('D','S')]
	tmp.df[match(grouse.TjD.list.mar[[gene]]$pop12$pos,tmp.pos),c('D12','S12')] <- grouse.TjD.list.mar[[gene]]$pop12[,c('D','S')]
	grouse.TjD.mar.unlist <- rbind(grouse.TjD.mar.unlist,tmp.df)
	} #ignore warnings

#calculate confidence intervals for D based on permutation on 1:S sites
D.CI.mat <- array(0,dim=c(4,10,20))
n <- 50
tmp.depth.loc <- grep('depth',names(grouse.mar.df3))

for(pop in 1:4){
	for(S in 1:dim(D.CI.mat)[2]){ #no of segregating sites
		for(rep in 1:dim(D.CI.mat)[3]){
			tmp.samp <- grouse.mar.df3[,tmp.depth.loc[pop]]>=n & (grouse.mar.df3[,tmp.depth.loc[pop]-1]/grouse.mar.df3[,tmp.depth.loc[pop]])>=0.01
			tmp.data <- grouse.mar.df3[tmp.samp,c(tmp.depth.loc[pop]-1,tmp.depth.loc[pop])][sample(1:sum(tmp.samp),S),]
			tmp.data$freq <- tmp.data[,1]/tmp.data[,2]
			dummy.count <- floor(n*tmp.data$freq)
			tmp <- (n-dummy.count)*dummy.count
			k.hat <- sum(tmp)/choose(n,2)
			a1 <- sum(1/seq(1,n-1))
			a2 <- sum(1/seq(1,n-1)^2)
			b1 <- (n+1)/(3*(n-1))
			b2 <- 2*(n^2 + n + 3)/(9*n*(n-1))
			c1 <- b1 - 1/a1
			c2 <- b2 - (n+2)/(a1*n)
			e1 <- c1/a1
			e2 <- c2/(a1^2 + a2)
			D.CI.mat[pop,S,rep] <- (k.hat - S/a1)/sqrt(e1*S + e2*S*(S-1))
			if(rep%%1000 == 0) cat('population ',pop,'; S ',S,'; rep ',rep)
			}
		}
	}
rm(pop,S,rep,tmp.data,dummy.count,tmp,k.hat,a1,a2,b2,c1,c2,e1,e2)
#run full version on euler
save(grouse.mar.df3,file='grouse_mar_df3.txt',ascii=TRUE)
#see tajima_bootstrap.R

#read in euler file
load('D_CI_mat.txt')

#some plots, 99% CI
plot(D12~S12,data=grouse.TjD.mar.unlist,pch=19,cex=0.3)
tst <- apply(D.CI.mat[4,,],1,function(X){quantile(X,0.995)})
lines(x=1:length(tst),y=tst,col='red')
tst <- apply(D.CI.mat[4,,],1,function(X){quantile(X,0.005)})
lines(x=1:length(tst),y=tst,col='red')

TD.boot.table <- array(dim=c(4,6,150))
for(i in 1:4){
	TD.boot.table[i,,] <- apply(D.CI.mat[i,,],1,function(X){quantile(X,c(0.0005,0.005,0.025,0.975,0.995,0.9995))})	}

grouse.TjD.mar.unlist2 <- grouse.TjD.mar.unlist
grouse.TjD.mar.unlist2[,c('pval1','pval2','pval3','pval12')] <- 1
tmp.pval.col <- grep('pval',names(grouse.TjD.mar.unlist2))
tmp.Dval.col <- grep('^D',names(grouse.TjD.mar.unlist2))
for(pop in 1:4){
	for(i in 1:nrow(grouse.TjD.mar.unlist2)){
		tmp.D <- grouse.TjD.mar.unlist2[i,tmp.Dval.col[pop]]
		if(!is.finite(tmp.D)){
			grouse.TjD.mar.unlist2[i,tmp.pval.col[pop]] <- NA
			next #jump out to next i
			}
		tmp.S <- grouse.TjD.mar.unlist2[i,tmp.Dval.col[pop]+1]
		grouse.TjD.mar.unlist2[i,tmp.pval.col[pop]] <- min(c(0.001,0.01,0.05,1,0.05,0.01,0.001)[c(min(c(1:4)[tmp.D < c(TD.boot.table[pop,1:3,tmp.S],100)]),max(c(4:7)[tmp.D > c(-100,TD.boot.table[pop,4:6,tmp.S])]))])

		}
	}

grouse.TjD.mar.summary <- aggregate(grouse.TjD.mar.unlist2[,-c(1,2)],by=list(grouse.TjD.mar.unlist2$cDNA),FUN=function(X){median(X,na.rm=TRUE)})
names(grouse.TjD.mar.summary)[1] <- 'cDNA'
grouse.TjD.mar.summary <- merge(summary.fst, grouse.TjD.mar.summary,all.x=TRUE)
#not a clear association between high fst and low D, though significant levels correlate a bit

fisher.test(table(grouse.TjD.mar.summary$pop.all.bootstrap.p<=0.05,grouse.TjD.mar.summary$pval12<=0.05))